Non-equilibrium electronic transport in a one-dimensional Mott insulator 
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We calculate the non-equilibrium electronic transport properties of a one-dimensional interacting 
chain at half filling, coupled to non-interacting leads. The interacting chain is initially in a Mott 
insulator state that is driven out of equilibrium by applying a strong bias voltage between the leads. 
For bias voltages above a certain threshold we observe the breakdown of the Mott insulator state 
and the establishment of a steady-state electronic current through the system. Based on extensive 
time-dependent density matrix renormalization group simulations, we show that this steady-state 
current always has the same functional dependence on voltage, independent of the microscopic 
details of the model and we relate the value of the threshold to the Lieb-Wu gap. We frame our 
results in terms of the Landau-Zener dielectric breakdown picture. Finally, we also discuss the 
real-time evolution of the current, and characterize the current-carrying state resulting from the 
breakdown of the Mott insulator by computing the double occupancy, the spin structure factor, and 
the entanglement entropy. 



I. INTRODUCTION 

The theoretical understanding of the non-equilibrium 
transport properties of strongly interacting systems in 
low dimensions has become a very active field of research, 
mainly due to the experimental activity in the fields of 
nanoscale materials^— and cold atomic gases,— as well 
as due to advances in theoretical methods designed to 
deal with both the non-equilibrium situation and elec- 
tronic correlations (see Refs. [H-0 for an overview and 
references therein). When considering non-equilibrium 
electronic transport, we have in mind a nanostructure 
that is subject to a large external voltage such that lin- 
ear response theory does not apply anymore. The main 
theoretical question that one would like to address is the 
dependence of the electrical current on the applied volt- 
age, i.e., the current- voltage characteristics, understand- 
ing not only the steady-state current reached on large 
time-scales, but also the transient regime appearing on 
shorter time-scales. Another important question is the 
characterization of the current-carrying state, contrast- 
ing its properties against equilibrium states in the ab- 
sence of a voltage. From the experimental point of view, 
knowledge of the full dependence of the electronic current 
on the bias voltage through an interacting nanostructure 
is a question of utmost importance, as this measurement 
is a standard technique to map out electronic energy lev- 
els and to observe many-body effects in nanostructures 
(see, e.g., Ref. HI for experimental work and Ref. on 
theoretical work). 

A paramount issue when studying transport in 



strongly interacting systems is the behavior of the in- 
sulating states characteristic of these systems, the most 
relevant of which is the Mott insulator (MI) state. Con- 
siderable theoretical efforts have so far been devoted 
to the study of non-equilibrium transport in nanos- 
tructures such as quantum dots (see, e.g., Refs. [Tol — 
l2lT ). Using state-of-the-art-numerical approaches, sub- 
stantial progress has been made in calculating the 
current- voltage characteristics and non-equilibrium prop- 
erties of some basic models, such as the interacting 
resonant level modelii or the single-impurity Anderson 
modeliii— — — — as well as in understanding their tran- 
sient behaviori 14 i 22 i 23 Whereas quantum dots with an 
odd number of electrons exhibit perfect conductance in 
the low bias regime due to the Kondo effect (21 an ex- 
tended region with repulsive interactions, an even num- 
ber of electrons, and at half filling is an insulator. The 
crossover from single quantum dots to this Mott insulat- 
ing state has been studied in Refs. I25W27I on the level of 
linear response theory, showing that the ground state al- 
ternates between a conducting state for an odd number of 
sites and an insulating state for an even number of sites. 
Of course, in the limit of large systems, the difference 
between N and N + 1 electrons becomes irrelevant and 
the perfect conductance in a system with an odd num- 
ber of electrons can only be observed at, with respect to 
experiments, unrcalistically low energy scales. 

In this paper, we shall thus turn our attention to non- 
equilibrium electronic transport through an extended in- 
teracting region, described by the one-dimensional Hub- 
bard model with repulsive onsite interactions. Specifi- 
cally we consider a one-dimensional (ID) system consist- 
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FIG. 1: (Color online) Sketch of the one-dimensional nanos- 
tructure described in the text, with an extended interacting 
region connected to non-interacting leads. Open (solid) sym- 
bols represent non-interacting (interacting) sites. The dashed 
line shows the voltage profile used to drive the system out of 
equilibrium. The voltage is homogeneous in the leads, and 
interpolates linearly between those values in the interacting 
region. 



ing of an interacting region of length L- mt connected to 
two non-interacting leads (see Fig. [IJ. The interacting 
region is initially in the MI state and we focus on the 
strongly interacting regime with interaction strength on 
the order or larger than the bandwidth. The sudden ap- 
plication of a large external voltage drives the system 
out of equilibrium and causes a time-dependent electri- 
cal current to flow through the interacting region, de- 
stroying the MI state. Our goals are first, to calculate 
the nonlinear current-voltage (I-V) characteristics; sec- 
ond, to contribute to the characterization of the current- 
carrying state; and third, to study the time-dependence 
of the entanglement entropy in this set-up. We employ 
the adaptive time-dependent density-matrix rcnormaliza- 
tion group (tDMRG) method ) 28 ' 29 which has been suc- 
cessfully used to compute the non-equilibrium dynamics 
of single quantum dotsi 11 ' 14 ' 23 ' 30 ' 31 

The problem of destroying a MI by subjecting it to a 
voltage or an electric field is currently attracting signifi- 
cant attention, both for the ll} 3 ^— and the 3D cases^ 
as well as in heterostructuresj 40 ' 41 In the case of an ex- 
tended interacting region, the voltage can be applied with 
different spatial profiles in theoretical simulations, giving 
rise to different physical mechanisms for destroying the 
MI state. In the set-up sketched in Fig. [TJ we linearly in- 
terpolate between the voltages set in the leads across the 
interacting region. As we shall argue, this gives rise to a 
many-body Landau-Zencr mechanism through which the 
MI breaks down. This picture has been advocated for 
in a series of studies by Oka and collaborators^ 2 - -34 ' 37 
who considered both a ring geometry pierced by a time- 
dependent flusS and a MI subject to a linear potential 
without including a coupling to \eadsM- Both approaches 
model the application of an electric field. One of their 
main results is that the breakdown of a MI is governed 
by the same physical laws as the one of a band insulator, 
with the difference that the band gap needs to be replaced 
by the charge gap of the strongly-interacting ML 32 ' 34 

Our setup is chosen to closely catch features of an ac- 
tual transport experiment by including the leads. Wc 
shall provide a qualitative comparison of our results with 



other cases recently addressed in the literature! 3 ' — 
Transport through extended interacting regions that are 
not necessarily in a MI state has been studied as well 
in Refs. 14214461 emphasizing as recurring themes the ap- 
pearance of nonlinear current-voltage characteristics and 
negative differential conductances. 

Our main result is the accurate numerical calculation 
of steady-state currents for the geometry of Fig. [TJ Wc 
find that the current-voltage characteristics can be de- 
scribed by an expression of the form 

J(V) =aVe- v ' /v , (1) 

in agreement with Refs. I32rl33l37l implying that at suffi- 
ciently large voltages, the system is driven to a conduct- 
ing state with J cx V. We show that V c cx A 2 , where A c 
is the Mott gap. In addition, we analyze several quanti- 
ties in the current-carrying state, with a particular focus 
on the double occupancy and spin correlations. While 
the current-carrying state still has a tendency towards 
antifcrromagnetic correlations, this instability is strongly 
suppressed compared to the MI state. However, neither 
the spin-structure factor nor the double occupancy, which 
is a measure of the interaction energy stored in the inter- 
acting region, saturate in the time window that we can 
access numerically. This suggests that the interacting 
region still undergoes a reorganization of internal energy 
while the particle flow in and out of the interacting region 
is already constant. The crossover from the insulating 
regime to the conducting regime is also reflected in the 
time-dependence of the entanglement entropy. We fur- 
ther show that this quantity behaves similarly to the case 
of global quenches: in our set-up, which is relevant for 
transport, the entanglement entropy increases linearly in 
time in the conducting regime. Here, the increase of en- 
tanglement is due to real particles moving around, differ- 
ent from the situation encountered in quantum quenches 
with homogeneous particle densities, in which propagat- 
ing collective excitations induce entanglement ! 47 ' 48 

One-dimensional Mott insulators can be realized ex- 
perimentally in several classes of materials. A promis- 
ing class of materials that have been suggested to realize 
ID MI are carbon nanotubes^ 9 - - — A recent experiment 
on carbon-nanotubc field-effect devices made from small- 
band-gap and nominally metallic carbon nanotubes has 
shown evidence for the realization of such a MI stated 
Theoretical wor k 49 ' 54 indicates that carbon nanotubes 
can be modeled by the Hubbard model on a two-leg lad- 
der geometry. Since in this effort we are interested in the 
generic behavior of a MI in the non-equilibrium regime, 
and since wc also need to keep the numerical effort at 
a manageable level, we will consider only ID chains, as 
opposed to ladders. Nevertheless our results may set the 
grounds for future studies on the appealing two-leg lad- 
der geometry. 

Besides realizations in nanostructures, the electronic 
properties of some quasi-one dimensional transition 
metal oxides are known to be well described by the one- 
dimensional Hubbard model. Most notably, Mott insu- 
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lator physics was found to be realized in SrCu02 and 
Sr 2 Cu0 3 and the specific question of the dielectric break- 
down of the MI state was experimentally addressed by 
Taguchi et al. in Ref. The actual physics of this ex- 
periment, however, may go beyond a simple Hubbard 
model description, as has been emphasized by Eckstein, 
Oka, and Werner^ 

An additional and related line of experimental research 
uses time-resolved photoelectron spectroscopy to drive 
systems with a gap into gapless phases (see, e.g., Ref. [57]). 
This method allows one to discriminate Mott insulators 
from other insulating states. 

The outline of the paper is the following. In section [TT1 
we present the model and briefly describe our numerical 
approach. In section IIIII we present our results for real- 
time currents, spin correlations, the double occupancy, 
and the entanglement entropy Section IIVI contains a 
summary and we discuss our results, contrasting them 
against the recent literature. 

II. MODEL AND METHODS 

To study the non-equilibrium transport in a Mott in- 
sulator we consider a one-dimensional chain with L elec- 
tronic sites. The chain is divided into three different re- 
gions: a non-interacting region at the left, representing a 
lead; an interacting region in the center, where the Mott 
insulator state is located; and another non-interacting re- 
gion at the right, representing another lead (see Fig. [TJ. 
This setup allows us to include the effects of the leads, 
complementary to the approach taken in Ref. |3~H . The 
number of sites of the left (right) lead is L\ (L t ), and in 
the interacting region is L[ nt . The Hamiltonian of the 
whole system can be written as 

H = H int + Hint-leads + Heads, (2) 

where 

H- mt = ~t" J2 (cLc z+la + h.c.) 

a,i—L\-\-l 

L\-\-Li nt Li+Li n t 

+e ^2 ^tr + U 2J n n n ii ( 3 ) 

(T,i— Li + 1 i — Li + 1 

is the Hamiltonian of a Hubbard chain with onsite 
Coulomb repulsion U > 0. t" is the hopping matrix ele- 
ment between the sites in the interacting region, and eo 
is the chemical potential in the interacting region. The 
second term in the Hamiltonian is 

Hint-loads = ~t''^2(c' LiCr CL l + la + h.C. 
a 

+c{ int+Lia c Lint+Ll+llT + h.c), (4) 



connecting the Hubbard chain to the leads with a hop- 
ping t' , resulting in a tunneling rate V = 2t' 2 . The third 
term in the Hamiltonian is 

Hlcads = -tlcads {c^d+ia + h.C.) 

cr,2— 1 

L-l 

-^loads {c\ a c i+ i a + h.c), (5) 

<J,i=L l +L illt + l 

where £i oa ds is the hopping matrix element in the leads. 
In most simulations, we set £' = t" and we use ii ea ds = 1 
as the unit of energy unless stated otherwise. In all the 
equations above c\ a represents the creation operator for 
an electron at site i and spin projection a ="f , -J-, rn a = 
4 a Ci a , and n.i = n lt + 714. 

We are interested in the time evolution of the MI state 
in the interacting portion of the chain when it is driven 
out of equilibrium by a strong voltage bias applied be- 
tween the leads. Therefore, we first need to find the 
ground state of the system when the interacting por- 
tion of the chain is at half-filling (eo = —U/2) and then 
solve the time-dependent Schrodinger equation for the 
perturbed system with this state as an initial condition. 
The former is accomplished by performing a ground-state 
DMRG^— calculation with N = L particles. To per- 
turb the system and to drive the chain out of equilibrium 
we add an extra term to the Hamiltonian Eq. ^j, which 
has the effect of adding an electric potential at time t = 0: 

L 

i=l 

where Q(t) is the Heaviside step function and 

f -V/2 i < L t 

Vi = { -(* - L c ) E for U < i < U + L int , (6) 
[ V/2 i>L t + Lint 

where L c = L{ + (L- ln t + l)/2. This mimics the effect of 
an electric field E = V/(L- m t + 1) acting in the interact- 
ing part of the chain, V being the bias voltage induced 
between the leads (see Fig. [1]). 

To solve the time-dependent Schrodinger equation we 
use the adaptive time-dependent DMRG techniqu o 28 ' 29 
with the methods introduced in Rcfs. [ll||ll2ll3al3l||6l| 
to simulate non-equilibrium transport. In some cases, 
we use systems with Li odd and L r even since we find 
that the finite-size effects in the currents are less severe 
for this configuration (compare with Refs. [23ll6lll62l for 
the case of few quantum dots). 

The tDMRG simulations are carried out using a third 
order Trotter-Suzuki breakup with a time-step of St = 
0.1/ti ea ds and under the constraint of a fixed, maximum 
discarded weight of 5p <~ 10~ 7 . In practice, this implies 
that one starts the time-evolution with a relatively small 
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number of states (to > 100) which then grows fast. The 
maximum number of states during the time-evolution is 
to = 1600 states. Since the accuracy of the numerical 
results solely depends on these control parameters, i.e., 
the discarded weight and the time step, tDMRG can be 
considered a quasi-exact method, as the numerical error 
can be estimated by varying St and 5 p. 

In non-equilibrium, the entanglement encoded in the 
time-dependent wave function is not bounded by any area 
law as is the ground-state entanglement^ and may in- 
deed increase extensively as a function of time. Typi- 
cally, in so-called global quenches (i.e., the instantaneous 
and homogeneous change of one parameter on all sites) 
one finds a linear increase of the entanglement entropy 
(the von- Neumann entropy) S v n ~ t with time (see, e.g., 
Ref. for the case of conformally invariant systems). 
Since the number of states to used in a DMRG calcula- 
tion scales as^S 



m oc e s " N , (7) 

reaching long time scales is an exponentially expensive 
computational task whenever S v n ~ t. Understand- 
ing the time-dependence of S v n itself in generic set- 
ups is thus an important objective to judge limitations 
and capabilities of tDMRG, besides the general and 
timely interest in its time-dependence in various kinds 
of quenches i 47 ' 48 ' 64 

The fact that the number of states increases monoton- 
ically with time defines a maximum time for each simula- 
tion as the time at which the number of states needed to 
keep the discarded weight under a fixed value bp exceeds 
the maximum of to = 1600. Then, for representative 
parameters, we perform several runs with different dp 
to assess and assure the numerical quality of the data, 
which ultimately determines the maximum time £ max at 
which the data for a given observable are still sufficiently 
reliable. 

We define the symmetrized tunnel current as the aver- 
age of the two local currents connecting the interacting 
region to the left and right leads: 



3 = yE'i/w 1 -^ 11 - - 

<7 

+ c L 1+ L int ,<T C £i+£in t +i,<r - h.c). (8) 

We will denote the time-dependent expectation value 
of the symmetrized current by J(t) = (j(t)) whereas 
the time-averaged current will simply be denoted as J. 
The currents are measured in units such that J/V = 2 
corresponds to perfect conductance, i.e., Go = 2e 2 /h 
(e = h = 1 in our work). Local currents (ji) on other 
bonds are defined accordingly. 



III. RESULTS 

The structure of this section is the following. First, 
we present the real-time data for the electric current and 
discuss the properties of the steady-state currents estab- 
lished after the dielectric breakdown of the Mott insula- 
tor takes place. Second, we analyze the current- voltage 
characteristics. As the main result of the paper we find 
a simple function to describe the current as a function of 
the bias voltage and the value of the Lieb-Wu gap asso- 
ciated to the initial Mott insulating state, similar to the 
results reported by Oka et a L 34 i 37 Third, we character- 
ize the current-carrying state in the interacting region by 
studying the time evolution of the charge and the current 
profiles, the double occupancy, and the spin-spin corre- 
lations. Finally, we discuss the time-dependence of the 
entanglement entropy. 

A. Real-time data and steady-state currents 

Figure [5] shows some examples of the real-time data 
for the symmetrized tunnel current obtained from our 
simulations for U/t" = 5 and two values of T. The tran- 
sient behavior, in general, can be expected to depend 
on both the tunneling rate, set by T = 2(t') 2 , and the 
voltage. For a small interacting region coupled to non- 
interacting leads, the transient regime has been studied 
in Refs. HmHjU-ill 

In our results, for all voltages, the generic behav- 
ior is that the current first goes through a transient 
regime, with a maximum reached in the time window 
< t < 1/r. The figure shows that the time-scale for 
reaching the first maximum is independent of the bias, 
while it clearly depends on F (this is obvious if one plots 
the results versus time in units of l/ii ea ds)- Then, accom- 
panied with oscillations whose period decreases with in- 
creasing voltage V, we reach a quasi-steady state regime 
(typically at times tT > 2. ..6) where the current is con- 
stant, apart from oscillations. The amplitude of the os- 
cillations decays as the steady-state is approached, yet 
from our data we cannot determine whether this decay 
is an exponential one or not. The period t Q of the oscil- 
lations is a monotonically decreasing function, similar to 
the case of single quantum dots in which t a oc 1/Ui 65 ' 66 

The time window over which the steady-state current 
can be sustained on a finite system can in principle de- 
pend on both L and Li nt . L trivially limits the accessible 
time-scales to t < t Tec = 2(L — Li ni )/vF, where vp is the 
Fermi velocity in the leads ; 14 ' 23 ' 61 since by that time, the 
perturbations induced in the leads by the application of 
the bias have traveled from the interacting region to the 
boundary and back, then perturbing the quasi-steady- 
state currents. Lj n t does not pose any limit on the sta- 
bility of the steady-state regime for the setup considered 
here, because the bias voltage is introduced locally as a 
homogeneous electric field. Therefore we choose the val- 
ues of L and L; n t to give a value of t rec similar to the 
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FIG. 2: (Color online) Current J(t) as a function of time for 
Lint = 20, U/t" = 5, t" = t' , L = 101 and (a) T = 0.08ii ead s 
(i.e. t' = 0.2t lcads ), (b) T = 0.32f lcads (i.e. t! = 0.4t leads ). 
In (a), Vyti ca ds = 0.5,1,1.4,2,2.5 (bottom to top) and in 
(b), V/t lcadB = 0.5, 1, 1.5, 2, 3, 4, 5 (bottom to top). Note that 
the maximum time reached in these simulations are (a) t = 
40/tie a da and (b) t — 20/ii oads , in units of the inverse hopping 
matrix elements in the leads. The arrows in (a) indicate the 
time interval used to compute the steady-state current J for 
V = 1.4i lcads . 



< max discussed in the previous section, t T 



B. I-V characteristics 

In this section we focus on the steady-state current and 
its dependence on the various parameters of the model, 
presenting results obtained from extensive numerical cal- 
culations. In practice, we compute the steady-state cur- 
rent by averaging over one or two periods of the oscilla- 
tions at the longest times reached in the simulations (but 
t < irec) t° reduce the effect of the oscillations. An ex- 
ample is shown with arrows in Fig.[2ja) for V — 1.4/i ca d s - 
We shall find that the current is a simple function of the 
bias voltage with all the microscopic details of the model 
encoded in the two coefficients a and V c in Eq. ([T]). I-V 
curves were previously presented for the ring geometry 
for very short chains and in that case, the currents were 
extracted from the short-time dynamics^ 

Figure [3] shows the steady-state current J as a func- 
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FIG. 3: (Color online) Current-voltage characteristics of the 
MI. Symbols represent the tDMRG results for the steady- 
state current, computed from the data shown in Fig. [2ja) as 
explained in Sec. lIIlAI Dashed curves are fits to the function 
J(V) = aVe~ v °/ v , a and V c being free parameters in the fit. 
(a) results for U/t" = 5 with t' = t" = 0.2ti cads . The plot 
includes data from two different system sizes to demonstrate 
that finite-size effects are small, (b) The inset shows the same 
data on a log-linear scale for L = 101. The agreement with 
the fit to Eq. |T|) is excellent (except for very low bias voltages, 
where J is of the order of our numerical accuracy) . 



tion of the bias voltage V for L- mt = 20 and U/t" = 5 
with t' = t" = 0.2£i oa d s . The data from our numerical 
simulations for J as a function of the bias voltage fit to 
Eq. (Q} with an excellent agreement, a and V c being the 
fitting parameters. Therefore, for values of V < V c below 
the threshold V c , J is exponentially suppressed whereas 
for values of V > V c above the threshold, J increases lin- 
early. The exponential term is dominant at low bias and 
causes the suppression of the current and represents the 
Landau-Zener tunneling rate^ across the Mott gap. The 
linear term is dominant at large bias and represents the 
motion of current-carrying excitations across the chain 
in the conducting regime. 

Figure [4] contains the TV curves for several different 
U/t", keeping t' and t" fixed. Motivated by Fig. 2 from 
Ref. [32|, we have plotted the steady-state current as a 
function of V/A^, where A c is the charge gap. We have 
calculated the charge gap for finite systems with £j nt = 
20 sites, not connected to any leads, using 

A c = [E (N + 2,S Z )+E (N- 2, S z ) - 2E Q (N, S z )}/2 , 

(9) 

where Eq(N, S z ) is the ground state energy in subspaccs 
with N fermions and a total spin projection S z . Using 
this, and by also plotting the current in units of U 2 , all 
curves collapse on a single one, which, in particular, sug- 
gests V c ex A 2 ., as expected for a Landau-Zener type of 
breakdown of the MI state. 32 ' 34 As we show here, this im- 
portant fingerprint of Landau-Zener physics also survives 
upon coupling the interacting region to leads. 
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FIG. 4: (Color online) I-V curves for several U/t" at fixed 
values of t' — ti ea ds and t" = 0.6fi ca ds- A c is the charge 
gap computed for an isolated chain of Lint = 20 from Eq. ([9]). 
Symbols are tDMRG data for L = 100 (L = 80 for U/t" = 7). 
Lines are guides to the eyes. 
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FIG. 5: (Color online) I-V curves for t' / t" at a fixed 

U/t" = 5 With t" = 0.6tlead S . f'Aloads = 0.4,0.6,0.8,1 

(bottom to top). Dashed curves are fits to the function 
J(V) = aVe~ Ve ^ v , with a and V c being free parameters in 
the fit. 



Wc here therefore find essentially the same dependence 
of V c on U as Oka et al.r2i24 namely V c oc A. 2 , but 
with incorporating the leads into the model. There are 
some differences, though. First, it should be noted that 
our time-averaged current is extracted from simulations 
that reach much longer times than Ref. |H where only 
the short-time dynamics was available to estimate the 
steady-state currents. Second, we do not find an abrupt 
increase of the current at the threshold voltage, in con- 
trast to Ref. Therefore, our data are in a better 
agreement with the result of mapping the problem to a 
quantum walk (see Fig. 3 in Ref. 1691 ). We attribute the 
quantitative differences between Fig. 5 in Ref. |32| and 
our Fig. [3] to the difference in the calculation of J, the 
fact that our systems are larger, and the inclusion of the 
leads. 

To further explore the effect of the leads on V c we have 
computed I-V curves for a fixed value of U /t" = 5 and 
several t' , as shown in Fig. [5] Qualitatively, a larger t' 
leads to an overall increase of the current as reflected 
in the i'-dependence of a to be discussed later on. The 
threshold exhibits a weak dependence on t' as well, as 
we demonstrate in Fig. [Ha). Our observation is that 
V c (t' < t") > V c (t' = t") and V c (t' > t") < V c (t' = t"). 
The latter behavior can be explained by the observation 
that close to the interface, the local charge gap depends 
on t': t' < t" leads to a slightly enhanced gap compared 
to the bulk gap and vice versa. As a consequence, the 
double occupancy (di) = (discussed in detail 

below) in the interacting region is enhanced close to the 
interface compared to the bulk value for t' > t", while 
it is suppressed for t' < t". Therefore, for t' < t" the 
contacts suppress the current, giving rise to an increase 
of V c . In the case of t' > t", the largest local gap is in 
the bulk of the MI and decreases towards the boundary. 



The decrease of V c as t' —> iicads can be understood as 
a consequence of a smaller mismatch between t',t" and 
Pleads m that limit, which should give rise to a increase in 
the transmission of electrons across the interface region. 
Note that we observe that boundary effects in the initial 
state typically decay to the bulk values over a distance of 
about 5 sites, suggesting that L- m t — 20 is a reasonable 
choice to probe both the bulk and contact properties. 

Next, we address the dependence of the prefactor a 
on t' . The coefficient a sets the value of the differential 
conductance in the conducting regime. We present our 
results for a and various combinations of t' in Fig. [H^b), 
in units of Go = 2e 2 /h. Interestingly, in all cases stud- 
ied, a < 2Go- Moreover, this coefficient a monotonically 
increases with t' or T = 2t' 2 . To summarize, a depends 
on both t 1 and U and, phcnomcnologically, we find that 
a oc U 2 results in a convincing collapse of the I-V curves 
for U > 4i" (compare Fig. [4}. 




0.4 0.6 0.8 1 0.4 0.6 0.8 1 
t'/t. , t'/t, . 

leads leads 

FIG. 6: (Color online) (a) Threshold voltage V c vs. t'\ (b) 
prefactor a in Eq. Q vs. t'. U/t" = 5, t" = 0.6t lcads , L = 100. 
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We have also studied the dependence of the I-V curves 
on L int (not shown in the figures). We find that 



V c oc (L int + 1)A;! and 



1/(L 



int 



i) 



(10) 



This suggests that the breakdown should be viewed as 
field-driven with E = V/ {L- lnt + 1) taking the role of the 
electric field. We may therefore rewrite Eq. (JlJ as: 



(11) 



J = aEexp(-E c /E) 



This interpretation is in agreement with Refs. I32ll34ll39l . 
and we stress that the functional form of the I-V curve 
described by Eq. (JTT|) holds despite the presence of the 
leads. As we have shown here, the effect of the leads is a 
small deviation of 5, and the threshold field E n from the 
bulk values (compare Fig. [HI and Refs. I32II3I 



C. Characterization of the current-carrying state 

The goal of this section is to characterize the current- 
carrying state in the interacting region. To this end, 
we measure the electronic density and electronic current 
density profiles in the interacting region, the average dou- 
ble occupancy, and also the spin-spin correlations, yield- 
ing the spin structure factor. 



1. Density and current profiles 

Figure [7^ a)- (c) show the charge density (n./) as a func- 
tion of position at different times for U/t" = 5, V = 
2iieads, an d L int = 20 and the corresponding local cur- 
rents (ji) in (d)-(f). In the steady state, the charge in the 
interacting portion of the chain has a linear profile fol- 
lowing the profile of the applied bias. The overall charge 
density in the Hubbard chain remains at half-filling. 

From the results for the local currents, we see that 
the currents take finite values on all sites, which actu- 
ally happens immediately after applying the potential. 
This clearly distinguishes the breakdown mechanism in- 
duced by a linear profile from other spatial forms of the 
bias voltage. For instance, in the simplest case in which 
Vi = in the interacting region and Vi = ±V/2 in the 
left (right) lead, the physics underlying the breakdown is 
quite different as we have verified in additional simula- 
tions (results not shown here). In this case, the redistri- 
bution of the charge inside the interacting region can be 
described as an effective doping of the MI region from the 
two interfaces. This implies that the bulk of the interact- 
ing region will experience the effects caused by turning 
on bias with a delay, set by the length of the interacting 
region. 

Turning back to Fig. [3 to justify that the steady state 
in an extended system has been reached, the currents 
need to be constant both in time and space. From Fig. [7J 
we see that (ji) = const is not fulfilled, although the 
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FIG. 7: (Color online) (a)-(c) Charge density (n/) as a 
function of position at different times tT = 0.8,1.6,2.4 for 
U/t" = 5, t' = t" = 0.2t loads , V = 2t lcads , and L = 100. 
(d)-(e) Profile of the local currents (ji) at the same times as 
in panels (a)-(c). 



charge flow in and out of the system is constant, apart 
from the relatively small oscillations discussed before 
(compare Fig. This suggests that for the time scales 
reached in our simulations, the interacting region still 
undergoes a reorganization of charges and local energies. 
Indeed, from the data of Fig. [Jj we find (ji — ji-i) / 0, 
even at tT ~ 2.5. 



2. Double occupancy 

Figure [5] shows the average double occupancy per site 
in the interacting portion of the chain 



d av (t) 



1 

Lint 



Ll +£int 

E (*(*)> 



(12) 



as a function of time for U/t" = 5 and V/t\ ca< i s = 
0.5,1,1.4,2. At all these voltages, the average double 
occupancy oscillates with a period given by t Q = t (V) 
that decreases with increasing voltage 1/T^, similar to the 
behavior of the currents. 

Depending on the bias voltage two different behaviors 
can be observed. For bias voltages below the threshold 
V < V c , i.e., in the regime of exponentially suppressed 
currents, d av (t) is essentially constant, apart from the 
oscillations. For bias voltages above the threshold V > 
Vc, d av (t) increases according to 



d av (t) = A + B t + C cos(D t) , 



(13) 



i.e., linearly in time after averaging over the period 
t Q = 2n/D. The slope B can be interpreted as the rate 
of the production of pairs of doublons and vacancies in- 
duced by the effective electric field i 34 ' 37 Quite notably, 



8 



the double occupancy never saturates over the time win- 
dow simulated, i.e., a steady-state regime for this quan- 
tity is not reached in our simulations, even if the system 
is in the steady-state regime for the tunneling current. A 
similar observation has been made in the DMFT study 
by Eckstein et aZ.,— who also report a monotonically in- 
creasing double occupancy d av (t) in the steady-current 
regime. They ascribe this to the fact that the work done 
by the field is proportional to j E, which in a regime of 
constant currents is a constant. Hence this increase in 
energy has to go into the internal energy of the MI, in 
the absence of any dissipation or leads. 

We shall here elaborate in more detail on this reason- 
ing, adopting it to our set-up that includes the leads. To 
explain the time-dependence of d av of Eq. (|13j) we exploit 
the fact that the equation of motion for the average dou- 
ble occupancy operator d av = (1/L int ) Yli=L,+T n it n 4 
is the same as the one for the interaction energy. After 
some straightforward algebra, one gets 



U/t"=5, t"=0.2t 



leads 



dt av [/L int 



-f E ■ 



(14) 



i=L L + l 



where T = Ti, lt + Tint-leads is the kinetic energy opera- 
tor involving sites at the interacting region, and E is the 
constant electric field. For times in the steady-current 
regime, the time integration of the second term in the 
RHS gives a linear dependence on time, as the current is 
approximately constant. Assuming that d av is small, as 
Fig. [7] suggests, we can expand the quantum mechanical 
average of the kinetic energy operator in the interacting 
region as (T int ) m T + edd av + 0(dl v ), where T is the 
kinetic energy of the filled lower Hubbard band, and e<i is 
the kinetic energy of a doublon. As a filled band cannot 
increase its kinetic energy, the time derivative approxi- 
mates as d(Ti nt )/dt e^dd^/dt. With this assumption 
one can move the contribution from T lnt to the LHS of 
Eq. (|14p and conclude that the time derivative of the 
average double occupancy is 



dt 



dt 



Tint — leads 

+ E :U + 0{d av {tf 



i=Li + l 



(15) 

where all operators have been substituted by their quan- 
tum mechanical averages, and we have changed the equal- 
ity in Eq. (| 15[) to a proportionality to accommodate the 
term stemming from the kinetic energy of the doublons. 
The first term in the RHS is the energy flowing out of 
the interacting region carried away by the particles trans- 
fered to the leads. If the interacting part is an isolated 
system as in Ref. l39l . this term is absent. The interpre- 
tation of Eq. (|15p is that although the establishment of 
the steady-current regime implies a linear increase of the 
double occupancy and therefore of the interaction energy, 
part of this energy is transfered to the leads when acceler- 
ated particles leave the interacting region. This reduces 
the rate at which the double occupancy increases, allow- 
ing the system to stay in the steady-current regime for a 




FIG. 8: (Color online) Average double occupancy d av (t) in 
the interacting region as a function of time for U/t" = 5, 
V/tioads = 0.5, 1, 1.4, 2, and t' = t" = 0.2ti cads . 



longer time. The increase in the double occupancy im- 
plies that the system is not in a true steady-state in the 
sense that there are observables that depend on time. 

As for the existence and the nature of a true steady- 
state, two scenarios are conceivable. Obviously, due to 
the bounded spectrum, the increase of the double oc- 
cupancy cannot go on forever, so eventually it has to 
saturate. An extreme case would be that d av takes its 
maximum value d max = 0.5 compatible with the sys- 
tem being at half filling on average. Consequently, the 
current would vanish in this case. Alternatively, the in- 
ternal energy could saturate at sonic time, reflected in 
dn V = const < 0.5 (where the RHS is the maximum pos- 
sible value assuming an average half filling of the inter- 
acting region). In that case, a finite current flow would be 
possible and the energy gain due to particles getting ac- 
celerated by the electric field would have to be balanced 
by an equal energy flow into the leads. In cither case, 
the reorganization of doublons may take longer than the 
time needed to reach the steady-state regime for the cur- 
rent. In particular, it is well known that the dynamics of 
doublons in one-dimensional systems with U > W where 
W is the bandwidth can be slow, if not even delayed by 
metastable regimes (see, e.g., Refs. I70U721 for ID systems 
and Refs. l73l for higher dimensions). This aspect has also 
been touched upon in Ref. l39l . 

Unfortunately, our simulations are restricted in the ac- 
cessible times, and we can thus not clarify this point, 
leaving it as an open question for future research. 



3. Spin-spin correlations 

The (longitudinal) spin structure factor can be com- 
puted from the spin-spin correlations by taking a Fourier 
transform e [L t + 1, L t + L int ]): 
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FIG. 9: (Color online) (a) Spin structure-factor in the inter- 
acting region as a function of time for U/t" = 5 (t' = t" = 
0.2fi oa< ja) and V — 2ti ca d s . (b) Sk=-n in the interacting region 
as a function of time for V/t\ ca ds = 0.5, 1, 1.4, 1.8, 2, 3 (top to 
bottom). 



s k = 



7^. 



-i(l — m)k 



(Si S m ) 



(16) 



l.m 



Figure EJa) shows the spin structure-factor at different 
times for U/t" = 5.0 and V = 2ti ca ds- The main feature 
is the survival of antifcrromagnetic correlations in the 
current-carrying state: the shape of the spin structure- 
factor remains qualitatively the same, yet the weight 
of the k = 7r instability decreases steadily with time. 
We therefore show Sk=n(t) for several bias values V in 
Fig. MJ°) as a function of time. The figure unveils that, 
similar to the case of the average double occupancy, a 
steady-state regime for this observable is not reached in 
our simulations, i.e., on the longest times reached and 
for the system sizes considered here. Similar to the lin- 
ear increase of the average double occupancy, the Sk=n (t) 
decreases linearly in time. 



D. Entanglement entropy 

The entanglement entropy is defined as 



S ,j 



N.. 



-trfp^ln^)] 



(17) 



where p x is the reduced density matrix of a block of the 
length x (counting from the left end of the chain). The 
reduced density matrix and its spectrum of eigenvalues 
is a key object in DMRG and the entanglement entropy 
is thus one of the easiest accessible quantities^ 



Let us begin by recalling some established analytical 
results on the entanglement growth in quantum quenches 
in systems with conformal invariance: In a global quench 
(i.e., the change of a parameter on all sites), S v n,x oc t 
(Ref. |47|) whereas in a local quench, S v n, x oc \n(t/to)' 48 ^ 4 
For the case of a global quench, this has been confirmed 
in numerous numerical calculations, mostly using DMRG 
(see, e.g., Refs.MM)- 

Our situation is different, since a parameter - the bias 
voltage - is changed on all sites, but with an explicit 
site dependence. Our results for S v n,x = S V N,x(t) are 
displayed in Fig. [TUJ Panel (a) shows S v n,x = SvN,x(t) 
vs. x for all possible cuts accessed in a DMRG run for a 
fixed value of V = 2ii ca ds at different times. The overall 
increase of S v n,x as a function of time is evident. 

The key question here is how the flow of particles in 
the conducting regime gives rise to an increased entan- 
glement between, say, the left lead and the rest of the 
system. In particular, we expect basically no increase 
in the insulating regime of bias voltages V < V c oc A^. 
To address this point, we plot S V N,x(t) with x = Li in 
Fig. [TOl b) for several bias voltages. Generally, we find 
that S v m. x = ct. The dependence of the prefactor c on 
bias voltage V is shown in the inset of Fig. HOT b): its de- 
pendence on V can be described by the same functional 
form as the tunnel current, namely: 



c oc V exp(-V c , vN /V) 



(18) 



In particular, we find that V c . v n ~ V c within the accuracy 
of our numerical simulations, where V c is the threshold 
voltage extracted from Fig. [3] This is consistent with the 
picture that entanglement is predominantly induced by 
propagating particles, in contrast to global quenches, in 
which (rij(f)) =const. 

While the observation of S v n,x oc t implies that the 
simulations carried out here become exponentially costly 
at long times, we note that in similar set-ups, namely 
the case in which a confining potential of a linear form 
is present in the initial state and its removal at t = 
is used to drive the time-evolution, a weaker logarithmic 
increase is found. Specifically, for the exactly solvable 
XX model, Eisler et al. report S v n,x oc ln(£)i2i Two 
main differences between their set-up and ours need to 
be pointed out. First, in our case, the application of the 
bias Vi destroys the MI state and drives the current flow. 
Conversely, in the set-up of Ref. H3, the initial state al- 
ready has an inhomogeneous particle density, implying 
that correlations in the initial state are already very dif- 
ferent from the respective ground state ones at the same 
filling. These open questions and observations call for a 
full analysis of the behavior of S v n.x in global quenches 
with site dependent changes of parameter, that we leave 
as a future project. 
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FIG. 10: (Color online) (a) Entanglement entropy S V N,x(t) 
vs. position x of the cut taken in the bipartition for times 
tT = 0,0.8, 1.6,2.4. (b) S vN ,x(t) vs - time t for x = Li= 41. 
This cuts the system across the left link that connects the left 
lead to the interaction region (V/ti ca d s = 0.5,1,1.4,2, from 
bottom to top ). Inset: slope c of S v N,x(t) = ct computed 
in the time-interval tT £ [0.25, 3]. The dashed line is a fit to 
c = a V exp(— Vc.vn/V). In both panels, U/t" = 5, t = t" = 
0.2ti ca ds ■ 



IV. 



DISCUSSION 



In this paper we have studied the dielectric breakdown 
of a Mott insulator state in a realistic model with an in- 
teracting chain connected to non-interacting leads. Our 
numerical results confirm that the steady-state current as 
a function of the applied voltage is, over a wide range of 
voltages, described by a simple universal function, with 
all the microscopic details of the model encoded in two 
coefficients related to the conductance in the metallic 
regime and the value of the threshold voltage. Our work 
further elucidates the influence of contacts to the leads 
on the I-V curve: the overall current is a monotonically 
increasing function in the inverse tunneling rate 1 /T and 
the threshold, on finite systems, also exhibits a weak de- 
pendence on the contacts. 

The dielectric breakdown of the one-dimensional Hub- 
bard model was studied under dissipative tunnelling into 
the environment introduced by a imaginary gauge po- 
tential in Ref. ffl\ . and upon the application of a strong 
electric field introduced by a gauge potential in a ring 
geometry in Refs. (32^3137 1. The main conclusion of 



the latter papers is that the dielectric breakdown of the 
Mott insulator can be understood in the same terms as 
the one in band insulators, with the only change that 
the band gap has to be substituted for with the Mott 
gap in the calculation of the Landau-Zener parameter 
(i.e., the threshold field). The time-averaged current in 
small Hubbard rings shows a collapse of the currents to a 
universal curve when the currents are plotted as a func- 
tion of the Landau-Zener parameter^ sharing the same 
qualitative traits as our Fig. [31 with a negligible current 
before the breakdown and a linear I-V characteristics at 
biases larger than the threshold. An important conclu- 
sion of our work is the confirmation that the mechanism 
of the dielectric breakdown corresponds to the Landau- 
Zener tunneling mechanism and this mechanism survives 
upon coupling the interacting region to leads. 

It should be noted that another very recent tDMRG 
study by Kirino and Ueda£& has adressed the destruc- 
tion of the MI state upon application of a strong voltage 
as well. There are important differences with our work, 
though. In Ref. HH, no leads are included, and the bias is 
applied as a step- function function to a homogeneous MI, 
measuring the local current on the central link. While the 
I-V curve also shows an activated behavior, it is not clear 
whether the MI is also destroyed through a Landau-Zener 
mechanism in the set-up of Kirino and Ueda. In partic- 
ular, they report V c oc A c , in contrast to the results by 
Oka et al. and ours (compare Fig. [4}. This illustrates the 
rich and various physical scenarios that can underlie the 
breakdown of an insulating state, depending on the way 
the bias is applied. 

We have also studied the conducting state that is 
reached after the breakdown. The spin-spin correlations 
remain antifcrromagnetic in the steady state. A decrease 
in the amplitude of the correlations is observed as the 
bias exceeds the threshold value. The conducting state 
can also be distinguished from the initial insulator by 
an increase in the double occupation. In other words, 
the electric field creates excitations as pairs of doublons 
and holons that can carry the current^ The production 
rate of these excitations should be reflected in the pro- 
duction rate of doubly occupied sites. Quite notably, the 
time-dependence of both the double occupancy and the 
spin-spin correlations implies that the interacting region 
is not in a true steady state yet, in which these quantities 
would become stationary as well. 

Finally, we have also computed the time-dependence 
of the entanglement entropy. This quantity increases lin- 
early with time in the conducting regime, implying that 
tDMRG simulations become exponentially expensive at 
long times. On the positive side, studying transport 
through single quantum dots or extended structures has 
qualitatively the same computational complexity, since 
in both cases, S v n jX oc t (unpublished results for one 
quantum dot, see Ref. Hil ). Therefore, going from sin- 
gle to many quantum dots is equally feasible with this 
method, in contrast to other state-of-the-art techniques 
such as time-dependent NRG^ or real-time QMCJ£ In 
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the former, the complexity scales with the dimension of 
the interacting region and in the latter approach, the dy- 
namical sign-problem is expected to become more severe 
for structures more complex than a single quantum dot. 
We have here demonstrated that tDMRG can success- 
fully be applied to compute I-V curves of extended sys- 
tems, complementing our earlier work on non-equilibrium 
transport in the single-impurity problem ! 14 i 30 ' 61 

While our numerical analysis of several properties of 
the current-carrying state should be helpful in better un- 
derstanding its properties, we acknowledge that a more 
intuitive picture of the non-equilibrium steady-state is 
still desirable. For instance, one would like to con- 
trast the current-carrying steady-state against effective 
ground-state reference systems, an approach which in 
certain non-equilibrium cases works quite welli^ More- 
over, the interesting concept of an effective temperature, 
often used in studies of quantum quenches with a relax- 
ation into a thermalized state (see Ref. [13 and references 
therein), should be further explored for current-carrying 
states. 



In conclusion, we have shown that the dielectric break- 
down of the Mott insulator can be understood in terms 
of the Landau-Zener mechanism using a realistic setup 
that matches the experiment since we include the leads. 
Furthermore we have been able to fully characterize the 
steady-state currents as a function of the bias voltage 
with a simple form, covering the whole range of voltages 
and microscopic parameters, that can be experimentally 
tested. 
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